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Ch Abstract 

a 

^ W e study the existence and stability properties of clusters of alternating charge vortices in Bose-Einstein condensates. It is illustrated 
O 'that such states emerge from cascades of symmetry-breaking bifurcations that can be analytically tracked near the linear limit of 
W ihe system via weakly nonlinear few-mode expansions. We present the resulting states that emerge near the first few eigenvalues of 
2 the linear limit, and illustrate how the nature of the bifurcations can be used to understand their stability. Rectilinear, polygonal 
Q and diagonal vortex clusters are only some of the obtained states while mixed states, consisting of dark solitons and vortex clusters, 
are identified as well. 
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1. Introduction 

Vortices are among the most striking characteristics of 
nonlinear field theories in higher-dimensional settings [1]. 
They constitute one of the remarkable features of superflu- 
ids, while playing also a key role in critical current densities 
and resistances of type-II superconductors through their 
transport properties, and are associated with quantum tur- 
bulence in superfluid helium [2]. Vortices appear also in a 
wide variety of fields, ranging from fluid dynamics [3] to 
atomic physics [4] and optical physics [5] . 

A pristine setting for the study of vortices at the 
mesoscale has emerged after the realization of atomic 
Bose-Einstein condensates (BECs). In this context, so- 
called matter-wave vortices were experimentally observed 
therein [6], by using a phase-imprinting method between 
two hyperfine spin states of a 87 Rb BEC [7]. This achieve- 
ment subsequently triggered extensive studies concerning 
vortex formation, dynamics and interactions. For instance, 
stirring the BECs [8] above a certain critical angular speed 
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[9-12] led to the production of few vortices [12], and even 
of robust vortex lattices [13]. Vortices can also be formed 
in experiments by means of other techniques, such as by 
dragging obstacles through the BEC [14] or by the nonlin- 
ear interference of different condensate fragments [15]. Not 
only unit-charged, but also higher-charged vortices were 
observed [16] and their dynamical instabilities have been 
analyzed. 

While a considerable volume of work has been dedicated 
to individual vortices and to vortex lattices, arguably, vor- 
tex clusters consisting of only a few vortices have attracted 
less interest. The latter theme has become a focal point re- 
cently, through the experiments involving two- vortex states 
(alias vortex dipoles) [17,18], as well as three- vortex states 
[19]. In Ref. [17], vortex dipoles were produced by drag- 
ging a localized light beam with appropriate speed through 
the BEC, while in Ref. [18] they were distilled through the 
Kibble-Zurek mechanism [20], previously proposed and re- 
alized for vortices in Ref. [21]. For the nonlinear dynamics 
of the vortices in the dipoles of Ref. [18], see also the very 
recent analysis of Ref. [22]. In Ref. [19], different types of 
three- vortex configurations were produced by applying an 
external quadrupolar magnetic field on the BEC. The prin- 
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cipal ones among them were an aligned vortex "tripole" 
with a vortex of one topological charge straddling two other 
oppositely charged vortices, and an equilateral triangle of 
three same charge vortices. On the theoretical side, few 
vortex states have been considered also in a number of 
works. It was shown, in particular, that vortex dipoles (con- 
sisting of a pair of vortices with opposite circulation) are 
fairly robust in BECs [23]. More elaborate states, such 
as dipoles, tripoles and quadrupoles, were considered in 
Refs. [24,25]. Dynamics of such few vortex states in the 
weakly-interacting limit were performed in Ref. [26] , while 
the recent work of Ref. [27] connected the vortex dipoles to 
the instability of dark soliton stripes; see also the impor- 
tant earlier work of Ref. [28]. 

In this work, we present a unifying analysis of the exis- 
tence and stability of vortex clusters (consisting of alternat- 
ing charge vortices), by corroborating theoretical investi- 
gations and numerical computations. Our study in Section 
2 will be based on the low-density limit of near-linear ex- 
citations, where we will illustrate how they emerge (bifur- 
cate through symmetry-breaking bifurcations) from states 
of the two-dimensional (2D) quantum harmonic oscillator. 
Then, such theoretically identified states will be continued 
via numerical computations in Section 3 to the strongly 
nonlinear regime. From these continuations, we will be able 
to infer numerous previously undiscovered vortex cluster 
states, and to elucidate their stability properties (as well 
as compare to the theoretical predictions). Lastly, in Sec- 
tion 4, we will summarize our findings and present some 
directions for future work. 



2. Model And Theoretical Analysis 

We consider a quasi-2D (alias "disk-shaped") conden- 
sate confined in a highly anisotropic trap with frequencies 
uo z and uo± along the transverse and in-plane directions, 
respectively. In the case uo± <C uo z and fi <C huo z (where 
ji is the chemical potential), and for sufficiently low tem- 
peratures, the in-plane part u(x,y,t) of the macroscopic 
BEC wave function obeys the following (2 + l)-dimensional 
Gross-Pitaevskii equation (GPE) (see, e.g., Ref. [4]): 



ihd t u = 



2m 



V(r)+g 2D \u\ 2 - n 



u, (1) 



where is the in-plane Laplacian, while the potential 
is given by V(r) = (1 / , 2)muo 2 _r 2 (where m is the atomic 
mass). The effective 2D nonlinearity strength is given by 
92D = gsD/V^TTCiz = 2\ph\ aa z hjO z , with g^D = 47rft 2 a/ra, 
a and a z = ^h/muo z denoting, respectively, the three- 
dimensional (3D) interaction strength, the s-wave scatter- 
ing length, and the transverse harmonic oscillator length. 
Equation (1) can be expressed in the following dimension- 
less form, 
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where the density |^| 2 , length, time and energy are respec- 
tively measured in units of [2\f2ji aa z )~ x , a z , uo~ x and Huj z . 
Finally, the harmonic potential is now given by V(r) = 
(l/2)Q?r 2 , with Q = uo±/uo z . From here on, all equations 
will be presented in dimensionless units for simplicity. 

Below, we will analyze the existence and linear stability 
of the nonlinear modes of Eq. (2). Notice that numerically 
the relevant nonlinear states will be identified as a function 
of the chemical potential \i by means of a fixed point (New- 
ton iteration) scheme over a rectangular two-dimensional 
grid with suitably small spacing. We will also explore the 
linear (spectral) stability of the obtained states by means 
of the Bogoliubov-de Gennes (BdG) analysis. The latter 
involves the derivation of the BdG equations, which stem 
from a linearization of the GPE (2) around the stationary 
solution uq(x, y) via the ansatz 



u = u (x,y) + a(x,y)e LUJ + b*(x,y)e 



j*t 



(3) 



where * denotes complex conjugate. The solution of the 
ensuing BdG eigenvalue problem yields the eigenfunctions 
{a(x, y), b(x, y)} and eigenfrequencies uo. Due to the Hamil- 
tonian nature of the system, if uo is an eigenfrequency of the 
Bogoliubov spectrum, so are — uo* and — uo*. Notice that 
a linearly stable configuration is tantamount to Im(cj) = 
0, i.e., all eigenfrequencies being real. It is important to 
mention that in what follows we only resolve the relevant 
eigenvalues up to 10 ~ 2 due to computational domain con- 
straints and therefore all eigenvalues smaller than 10 ~ 2 will 
be omitted in the figures. 

Within the BdG analysis, a relevant quantity to consider 
is the norm x energy product of a normal mode with eigen- 
frequency cj, namely, 



E 



dxdy(\a\ 2 - \b\ 2 )uo. 



(4) 
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The sign of this quantity, known as Krein sign [29] , is a topo- 
logical property of each eigenmode. In particular, if this sign 
is negative and such a mode becomes resonant with a mode 
with positive Krein signature then, typically, complex fre- 
quencies appear in the excitation spectrum, i.e., a dynam- 
ical instability arises [29] . We refer to this as an oscillatory 
instability. Furthermore, dynamical instabilities may arise 
due to a real mode eigenfrequency becoming imaginary. 
This typically coincides with a bifurcation of a new state. 

In the context of Eq. (2), it is useful to consider the low- 
density (linear) limit. There, eigenstates of the 2D quan- 
tum harmonic oscillator arise in the form u nm (x,y,t) = 
ex.p(—iiJ J t)H n (x)H m (y) (as well as linear combinations 
thereof), where fi = Q(n + m + 1) and n, m quantify 
the order (and number of nodal lines) in each direction. 
This produces a linear limit whose first excited state has 
li = 20, and linear eigenstates u\q and uq\. Notice that one 
of their interesting linear combinations is u\q + ^oi? which 
creates the single-charged vortex even at this linear limit; 
this state exists for all higher values of \i and is shown in 
Fig. 1, but we will not be concerned with it further herein, 
as our focus will be on clusters of vortices. 
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Fig. 1. (Color online) Top left panel (a): Number of atoms as a function of the chemical potential for the different states for Q = 0.2 bifurcating 
from the dark soliton stripe. Bottom left panel (b): Corresponding atom number difference with respect to the dark soliton stripe branch. 
Bifurcation of vortex multipole states arises when their atom number difference from the dark soliton stripe vanishes. The corresponding 
theoretical predictions of such bifurcations are shown by the vertical lines. Right panels (from left to right and top to bottom): Density and 
phase profiles (left and right subpanels respectively) corresponding to bifurcating states from the dark soliton stripe (bottom right panel): 
the single vortex state, the aligned two- (vd), three- (3v), four- (4v) and five- vortex (5v) states. 



Each of the above mentioned linear eigenstates, uio and 
uqi, represents a "stripe" i.e., a state with a nodal line. 
These can be continued for higher chemical potentials \i. As 
\i increases, this state develops into a one-dimensional (ID) 
dark soliton stripe, which is an exact analytical solution of 
Eq. (2) in the absence of the trap [30]. However, it is well- 
known that such a state is dynamically unstable towards 
decay into vortex structures [31-33]. This decay can be 
understood from a symmetry-breaking bifurcation point of 
view [27,34]. In particular, it is possible to consider a two- 
mode (Galerkin-type) expansion, similar to the one used in 
the literature of double- well potentials (see e.g. Ref. [35]) 
in the form: 

u(x,y,t) = co(t)<f> (x,y) + ci(t)0i(x, y), (5) 

where co(t), c\(t) are complex time-dependent prefactors, 
while 0o = v>io(x,y), 4>i{x,y) = wo m and m > 1. The re- 
sulting equations and analysis are formally equivalent to 
the ones derived in Ref. [35] (see Eqs. (4)- (5) therein), with 
appropriate modifications of the inner products, but also 
with a fundamental difference. In the ID double- well set- 
ting, only symmetry-breaking bifurcations of asymmetric 
real solutions are predicted (the so-called 7r-states that have 
recently been experimentally observed in Ref. [36]). The 
richer 2D case enables bifurcations even when the relative 
phase Acj) between the complex order parameters Co and c\ 
is 7r/2. In particular, such bifurcations are generically pre- 
dicted at an atom number: 



where cjo, are the linear state eigenvalues corresponding 
to and 0i, while Iq = J (\>\(\)\dxdy and 1\ = J (\)\dxdy\ 
the critical chemical potential is given by /i cr = ujq + IiN cr . 



3. Numerical Results and Comparison with Theory 

The above mentioned two-mode theory provides explicit 
predictions for the bifurcation not only of the vortex dipole 
(vd) state when m = 2, but also for an aligned vortex tripole 
(3v) state (cf. the experimental observations of Ref. [19]) 
for m = 3, for an aligned vortex quadrupole (4v) state 
with m = 4, etc. In fact, there is an entire cascade of such 
bifurcations, as m increases, which occur progressively at 
/i^ r d = 10ft/3, /4 V = 860/19, a4 v = 890^/157, and ^ = 
726^/107 for m = 2,3,4,5 etc., respectively. Notice that 
the number of vortices of the resulting cluster is evident by 
the number of intersections of the single nodal line of uio 
with the m perpendicular nodal lines of wo m , as well as the 
7r/2 relative phase of their complex prefactors at these m 
intersections. Also, it is evident that that the sign chang- 
ing of uom at these intersections leads to an alternation of 
the ensuing vortex charges. Importantly, general bifurca- 
tion theory can be used to identify the stability character- 
istics of the resulting states. In particular, since the stripe 
is dynamically stable as it emerges from the linear limit, 
the vortex dipole state that arises from it upon the first 
symmetry-breaking "event" (m = 2) should inherit this 
stability. However, now, once the stripe has become unsta- 
ble, all higher bifurcations with m > 3 will necessarily re- 
sult into dynamically unstable states. 

Numerical results on the symmetry-breaking bifurca- 
tions resulting in the emergence of vortex cluster states 
from the first excited state (single dark soliton stripe) 
are summarized in Fig. 1 for Q = 0.2. The emergence of 
1 x m (m = 2,3,4,5,...) states can be observed to oc- 
cur respectively at (0.68,0.98,1.26,1.54) while the corre- 
sponding theoretical predictions are (/i^, A^r > Mcr ? Mcr) = 
(2/3, 0.91, 1.13, 1.36). Clearly, the two- mode approach cap- 
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Fig. 2. (Color online) Top left panel (a): Number of atoms as a function of the chemical potential for the different states bifurcating from 
the two dark soliton stripes (black lines) and from the dark soliton cross (orange [gray in printed version] lines) for fl = 0.2. Bottom left 
panel (b): Corresponding atom number difference with respect to the two dark soliton stripes branch. The occurrence of bifurcations at zero 
crossings of AN and the corresponding theoretical vertical line predictions thereof are similar to the previous (and following) figures. Middle 
column of panels (from top to bottom): Density and phase profiles (left and right subpanels respectively) corresponding to bifurcating states 
from the two dark soliton stripes (top): the aligned six- (6v2) and eight- vortex state (8v2). Right column of panels (from top to bottom): 
Density and phase profiles (left and right subpanels respectively) corresponding to bifurcating states from the X-shaped dark soliton cross 
(top): the diagonal six- vortex state (6x), the four and doubly-charged vortex state (5x), and the diagonal eight- vortex state (8x). 



tures the fundamental phenomenology, although a slight 
progressive degradation of the agreement on the critical 
point arises due to the increasing departure from the lin- 
ear limit. Additionally, the linear stability properties of 
the resulting states directly reflect the theoretical stability 
expectations discussed above; see also Ref. [34]. 

Importantly, this approach is not restricted to the first 
excited state. The advantage of the bifurcation method and 
of the wealth of states that can be derived from it is un- 
veiled, e.g., when considering the next set of excited states, 
namely the combinations with n + m = 2, i.e., the degener- 
ate states 1^20, U02 and u\\. In this setting, already a vortex 
quadrupole can be formed at the linear limit as U20 + 11102 
[37]. Interestingly, so can a doubly-charged vortex through 
U2o — uo2 + 2iun. However, these states do not present sym- 
metry breaking bifurcations and, therefore, are not consid- 
ered further in what follows. Focusing on the states that 
do, some prototypical examples are (i) the solitonic state 
consisting of two stripes (i.e., a two-dark-soliton state), i.e., 
the U20 state, (ii) an X-shaped dark soliton cross emerg- 
ing from U20 — ^02, as well as (iii) a ring dark soliton state 
[38] (see also Refs. [39,40]) arising from U20 + ^02- The U20 
state is one for which the generalization of the phenomenol- 
ogy of Fig. 1 is most straightforward as with 0o = ^20 and 
4>i = uom states with 2 x m (2 lines of m vortices each) are 
formed with m > 2. For example, for <pi = U03 and <j>i = 
1^04, such bifurcations are predicted at /i^J 2 = 2830/67, and 
/ij?J 2 = 29650/551 (see thick black vertical dashed lines in 
Fig. 2), respectively, leading to 6 = 2 x 3 and 8 = 2x4 



two-line vortex clusters. On the other hand, the X-shaped 
dark soliton pair is subject to similar symmetry-breaking 
bifurcations/destabilizations due to <pi = ^03, ^13, ^04 
etc. These bifurcations, occurring in turn at = 330/7, 
f4? = 990/17, and ^* = 3690/59 (see thin orange [gray 
in printed version] vertical dashed lines in Fig. 2), lead to a 
diagonal state with 6 vortices, a doubly-charged vortex in 
the center together with a four- vortex quadrupole around 
it, and an 8- vortex cluster of near-diagonal vortices. Finally, 
the ring dark soliton gets mixed with states of the form = 

qi(r) sm(W), where qi(r) = \J exp(— Or 2 /2), again 
for I > 3. Remarkably, these symmetry-breaking events, 
which can be predicted to occur, e.g., at /x^J = 50 and 
fj^l = 590/9 (see, respectively, orange [gray in printed ver- 
sion] and blue [dark in printed version] vertical dashed lines 
in Fig. 3), for I = 3 and 4, give rise to polygonal vortex con- 
figurations with the vortices now placed on the periphery of 
the circle. This way, vortex hexagons, octagons, decagons, 
etc. can be systematically constructed at will. 

In Figs. 2 and 3, we explore numerically the above 
symmetry-breaking events for the second excited states 
and show besides the bifurcation diagram also typical 
profiles of the waveforms, as they result from varying the 
chemical potential. Specifically, Fig. 2 depicts the bifurca- 
tions and profiles for the states arising from the two dark 
soliton stripes (middle column of panels) and the dark 
soliton cross (right column of panels). From the two-stripe 
soliton bifurcates a six-vortex (6v2) (2 x 3) state and an 



4 




0.6 0.8 1 1.2 1.4 1.6 1.8 



Fig. 3. (Color online) Top left panel (a): Number of atoms as a function of the chemical potential for the different states bifurcating from 
the ring dark soliton for fl = 0.2. Bottom left panel (b): Corresponding atom number difference with respect to the ring dark soliton. Right 
panels (from left to right and top to bottom): Density and phase profiles (left and right subpanels respectively) corresponding to bifurcating 
states from the ring dark soliton state (top left): the hexagonal (6r) and octagonal (8r) vortex states. Bottom right panel: vortex quadrupole 
state (note that this state does not bifurcate from the ring dark soliton). 




Fig. 6. (Color online) Top left panel (a): Number of atoms as a function of the chemical potential for the different states bifurcating from 
H = 4f2 for Q = 0.2. Bottom right panel (b): Corresponding atom difference with respect to the three-soliton state. Right panels (from left to 
right and top to bottom): Density and phase profiles of the nine- vortex state, the vortex-soliton-ring, the three-soliton state and the mixed 
state consisting of a dark stripe and six vortices. 



eight-vortex (8v2) state (2 x 4) in the considered regime. 
From the X-shaped dark soliton cross bifurcates a diago- 
nal six- (6x) and eight- vortex (8x) configurations, as well 
as a state with a vortex quadrupole surrounding a central 
doubly-charged vortex (5x). We also mention in passing 
that interesting additional bifurcation events (collisions 
and disappearances into "blue sky" bifurcations) arise, e.g., 
between the six- vortex cluster deriving from the X-shaped 
dark soliton cross and that deriving from the two-stripe 
soliton (see blue circle denoted by A in Fig. 2). On the other 
hand, in Fig. 3 we depict the states that bifurcate from 
the vortex ring: vortex hexagons (6r) and octagons (8r). 
In this figure we also depict the vortex quadrupole state 
(bottom right panel) which does not bifurcate from the 



soliton ring. In all the cases, good quantitative agreement 
is found on the critical points predicted by the theory and 
those observed numerically. 

Figures 4 and 5 complement the above existence picture 
with a systematic linear stability analysis of each of the cor- 
responding states for Figs. 2 and 3 respectively. The two- 
stripe soliton is dynamically stable near the linear limit, a 
stability inherited by its first symmetry-breaking offspring, 
the six- vortex (2x3) state; while the aligned eight- vortex 
case (2 x 4) is unstable from its existence onset (see left col- 
umn of panels in Fig. 4). Similarly, the X-shaped dark soli- 
ton cross state also bears imaginary eigenfrequencies for all 
values of \i and hence its derivative states inherit its dynam- 
ical instability (see right column of panels in Fig. 4). On the 
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Fig. 4. (Color online) BdG spectra for the two-dark-soliton stripe 
and its bifurcating states (left column) and for the dark soliton cross 
and its bifurcating states (right column) corresponding to the states 
described in Fig. 2. The inset in each of the spectra depicts a typical 
configuration at the chemical potential indicated. In all BdG spectra 
we depict the real (top subpanel) and imaginary (bottom subpanel) 
parts of the eigenfrequencies as a function of the chemical potential fi 
using the following color coding for the online (respectively, in print) 
version: blue (black) ^positive energy (Krein sign) modes, orange 
(dark gray)=negative energy modes, green (gray)=oscillatory unsta- 
ble modes, and pink (light gray)=non-oscillatory unstable modes. 
The presence of non-vanishing imaginary parts is an indication of 
instability. Note that the only stable solutions are the two-dark soli- 
ton stripe (top left panel) and its first symmetry breaking offspring, 
the six- vortex state (left panel of second row), which are stable only 
for low enough chemical potentials. 

other hand, the ring dark soliton is, as was also found ear- 
lier [37-40], always unstable, hence the polygonal vortices 
that derive from it inherit this instability (see Fig. 5). How- 
ever, it should be noted that the instability of such states 
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Fig. 5. (Color online) BdG spectra for the states bifurcating from the 
vortex ring and for the vortex quadrupole (bottom right) that does 
not bifurcate from the vortex ring corresponding to the bifurcations 
depicted in Fig. 3. For an explanation of the color codes see Fig. 4. 
From the above it is clear that the most robust state is the vortex 
quadrupole (which is dynamically stable for small \i except for a 
narrow interval of oscillatory instabilities). 
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Fig. 7. (Color online) BdG spectra for the bifurcating states depicted 
in Fig. 6. From left to right and top to bottom: the nine- vortex state, 
the ring dark soliton plus vortex, the three-dark-soliton stripe (which 
is stable for low values of the chemical potential), and the mixed 
soliton-vortex state. For an explanation of the color codes see Fig. 4. 

weakens as the chemical potential \i increases. Note that, 
the vortex quadrupole is stable in the linear limit. Since no 
state bifurcates from the vortex quadrupole this stability 
is generally preserved (apart from a small window of oscil- 
latory instabilities). 

A general comment on the observed BdG spectra is due 
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here. We can see that the spectra bear a large number of 
positive Krein sign modes (see blue (black) points is all 
BdG spectra) which appear to be asymptoting to appro- 
priate corresponding values in the large chemical potential 
limit. In addition, there is a number of negative Krein sign 
eigenfrequencies (see orange (dark gray) points is all BdG 
spectra) which may lead to oscillatory (upon collision with 
positive sign ones) or exponential instabilities. Generally, 
we can comment that the positive Krein sign eigenmodes 
correspond to the ground state "background" on which a 
particular solitonic or vortex (or mixed) state exists. Fur- 
thermore, these eigenfrequencies have a well-defined limit 
when fi is large as discussed e.g. in Ref. [41] (see also ref- 
erences therein). On the other hand, the negative Krein 
sign modes reflect the excited state nature of the consid- 
ered soliton or vortex (or mixed) states and are the ones 
which bear the potential for dynamical instabilities. While 
for solitonic states, we do not have a precise count on the 
number of the latter eigenmodes, in the case of multi- vortex 
clusters consisting of individual vortices an upper bound 
on the maximal order of such (negative Krein) eigenstates 
can be given by the number of vortices in the configuration. 
This, in turn, also gives an upper bound on the number of 
potentially observable unstable eigenmodes. 

Importantly, the present approach can be generalized to 
higher excited states. As merely a small sample of further 
exotic configurations that can emerge (some of which can 
even be structurally robust), we mention a 9- vortex cluster 
(a square grid of 3 x 3 vortex "particles" ) , which emerges 
from the linear limit as uso + iuos and can, therefore, be 
regarded as a higher excited analog to the single vortex 
and the vortex quadrupole states and is expected to be 
stable in the vicinity of that limit (at least with respect 
to purely imaginary [i.e., non-oscillatory] instabilities, see 
below). On the other hand, there are bifurcations from that 
state including the bifurcation of a ring with a vortex in 
its center (which was again discussed e.g., in Ref. [40]). 
Another state that exists and should be robust near the 
linear limit is a three-soliton-stripe. In fact, from such a 
state bifurcations again emerge due to the mixture with 
states such as 4>\ = UQ rn with m > 4, but also with 0i = 
uim with m > 3. We focus briefly on the latter, which is 
theoretically predicted to occur at /i^J s = 1250^/239 (see 
vertical dashed line in Fig. 6), because it gives rise to yet 
another novel type of state, namely a mixed state between 
vortices and a dark soliton: this soliton- vortex mixed state 
has a dark soliton stripe nodal line plus two additional lines 
over each of which three vortices reside. Similar states with 
8-, 10- etc. vortices beside the soliton also exist, arising 
through subsequent bifurcations. 

To corroborate the above theoretical predictions we show 
in Figs. 6 and 7 (again for Vt = 0.2) some prototypical ex- 
amples of states that emerge from the third excited linear 
branches and their corresponding BdG stability spectra. 
The 9- vortex "crystal" is stable close to the linear limit in 
the sense that its BdG spectrum possesses no imaginary 
mode. However there is a complex mode inducing an os- 



cillatory instability. At fi w 0.95 the vortex plus ring soli- 
ton state bifurcates from the 9- vortex "crystal" inducing a 
complex mode. However, the latter is small and vanishes 
again at \i « 1.6. The BdG spectrum of the vortex plus 
ring soliton state contains no purely imaginary mode but 
many complex modes. Note that the eigenvalue spectrum 
looks fairly similar to the ring dark soliton case, but the 
imaginary part of the modes is not due to imaginary modes 
which are created by bifurcations. In this case, modes with 
positive energy cross zero and thus obtain a negative en- 
ergy. These negative energy modes then collide with modes 
with positive energy and create the complex modes. More- 
over, we show the three-soliton-stripe, as well as the mixed 
soliton-six-vortex cluster state emerging from its first bi- 
furcation event in Fig. 6. The three-soliton state is stable 
near the linear limit (of /x = 0.8) but becomes destabilized 
at fi = 1.05 (and then further so at 1.06, 1.1 and 1.3). 
The first bifurcation gives rise to the very weakly unsta- 
ble soliton- vortex state predicted above; the critical point 
for this bifurcation is found to be at \i = 1.034, once again 
in very good agreement with the full numerical result. The 
present approach can naturally be extended to a multitude 
of additional states which, however, are expected to be dy- 
namically unstable; therefore, having presented the most 
fundamental ones, we will not proceed further with such 
considerations. 

4. Conclusions and Future Directions 

In conclusion, in the present work we have shown that a 
detailed understanding of emergent vortex cluster states is 
possible through a near-linear-limit approach. This involves 
identifying the possible linear states and tracking system- 
atically the symmetry-breaking bifurcations that can arise 
from them. This allowed us not only to discuss a cascade 
of bifurcations from the first excited state (in the form of 
aligned vortex clusters), but to also reveal a broad class of 
states emerging from the second- and third-excited states. 
These were not only rectilinear states (with "soliton-type" 
stripes), but also soliton rings and rings with vortices, vor- 
tex quadrupoles or nine- vortex-crystals, as well as various 
clusters of vortices that derive from some of these states, in- 
cluding vortex hexagons, octagons, decagons, nx m states 
(of m vortices sitting along n stripes), soliton- vortex states, 
and so on. 

We were also able, based on the general bifurcation struc- 
ture of the problem, to reveal which ones among these states 
are expected to be most robust, such as the vortex dipole 
or quadrupole, and some of the emerging vortex clusters, 
such as the 2 x 3 or the soliton- vortex state. 

It would be especially interesting to extend this picture 
to different dimensions. On the one hand, one can con- 
sider effects of anisotropy (by changing the strengths of 
the two different trapping directions). This should enable a 
"dimension- transcending" picture, as extreme anisotropics 
should allow to observe how the system transitions between 
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quasi-one- and genuinely two-dimensional dynamics. On 
the other hand, it would be particularly relevant and inter- 
esting to attempt to extend such considerations into three- 
dimensional settings, and understand how relevant ideas 
generalize and potentially give rise to structures such as 
vortex rings that have been observed in pertinent experi- 
ments [33,42]. Such investigations are currently in progress. 
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